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Two-dimensional hydrodynamic simulation of an accretion 
flow with radiative cooling in a close binary system 
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ABSTRACT 

Two-dimensional numerical simulations of an accretion flow in a close binary system 
are performed by solving the Euler equations with radiative transfer. In the present 
study, the specific heat ratio is assumed to be constant while radiative cooling effect 
is included as a non-adiabatic process. The cooling effect of the disc is considered by 
discharging energy in the vertical directions from the top and bottom surfaces of the 
disc. We use the flux-limited diffusion approximation to calculate the radiative heat 
flux values. Our calculations show that a disc structure appears and the spiral shocks 
are formed on the disc. These features are similar to that observed in the case of an 
adiabatic gas with a lower specific heat ratio, 7 = 1.01. It is found that when radiative 
cooling effect is accounted for, the mass of the disc becomes larger than that assuming 
7 = 5/3, and smaller than that assuming 7 = 1.01. It is concluded that employing 
an adiabatic gas with a lower specific heat ratio is almost a valid assumption for 
simulating accretion disc with radiative cooling effect. 

Key words: accretion, accretion disc - hydrodynamics - methods: numerical - bina- 
ries: close - radiative transfer 



> . 1 INTRODUCTION 



The angular momentum loss of the gas in the accretion disc 
in a close binary system is one of the very important phe- 
nomena in the astrophysics. The most accepted model that 
explains the mech anism of angular momentum transport i s 
the a disc model (IShakura fc SunvaevllT97ol |Pringldll98lD . 
In this model, viscosity originating from turbulence, mag- 
netism or whatever, is supposed to transport the angular 
momentum. However, in this model the details of how the 
viscosity is generated are still unknown. 

An alternative model that needs no viscosity for an- 
gular momentum transfer is the spiral shock model in 
which spiral shaped shock waves are formed in the accre- 
tion disc due to the tidal f orce of the mass losing star. 
This m odel was proposed by ISawada. Matsuda fc Hachisul 
lll986alh They conducted two-dimensional hydrodynamic 
calculations of a flow of invisicid adiabatic gas, and found 
the spiral shocks numerically. Since then a number of 
two-dimensional simulations have been carried out, and 
it has been claimed that t he spiral shocks really ap- 
pear in the accretion discs llSawada. Matsuda fc Ha chisu 
Il986bl: ISawada et alJll987t ISpruit et alJll987t lLarsonlll988t 



* E-mail:junichi@cfd. mech.tohoku.ac.jp (JS) 
f E-mail:sawada@cfd. mech.tohoku.ac.jp (KS) 
\ E-maihohnishiOcfd. mech.tohoku.ac.jp (NO) 



| SpruitJll98sl: iRozvczka fc Spruidll98ll iMatsuda et alJ|l990( : 
Savoniie, Papaloizou & Li nl 11994) . Besides. ISpruitJ dl987D 
obtained self-similar solutions of spiral shocks for idealized 
accretion discs. 

For this spiral shock model, two major points have been 
questioned as regard to the validity of the model: 

(i) The three-dimensional effe ct might suppress the forma- 
tion of spiral shocks Il989t iLin. Papaloizou fc Savoniie] 
ll990dl3lL~ubow fc Pringlelll993lh 

(ii) The cooling effect due to radiation is not taken into 
account explicitly. Though such effect has been considered 
by assuming a lower specific heat ratio for adiabatic gas, the 
obtained disc was still too hot. 

As regard to the question (i) , three-dimensional simula- 
tions have been carried out mostly by particle methods. Such 
simulations tend to give rather diffusive results if a sufficient 
number of particles are not available. Spiral shocks are likely 
disappeared particularly for those cases with a lower spe- 
cific heat ratio for which the pitch angle of the spira l shock 
decreases. For example, lYukawa. Boffin fc Matsudal (Il997t) 
carried out three-dimensional simulations using smoothed 
particle hydrodynamics method. As many as (5 — 6) x 10 4 
particles were employed in the calculations. They showed 
that an accretion disc was formed and spiral shocks appeared 
for the case of 7 = 1.2, but spiral shocks disappeared for the 
case of 7 = 1.1 and 7 = 1.01. 
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Three-dimensional simulation using either the finite 
dif ference or finite volume m ethod was first carried out 
by Sawad a fc Matsudal (^292;). They calculated the case of 
7 = 1.2 with a mass ratio of unity using the upwind TVD 
Roe scheme and a generalized curvilinear coordinate sys- 
tem. Although the calculation was continued up to only a 
half revolution period because of the limitation of CPU time, 
they successfully showed the existence of spiral shocks on the 
accr etion disc. 

iBisikalo et al l Jl997ah. iBisikalo et all |l997bft . 



[Bisikalo et alJ <1998aft . IBisikalo etafl (1199811 " and 
Bisik alo et alJ il998a) carried out three-dimensional 
numerical simulations of accretion discs by means of the 
finite difference method. They employed a TVD Roe scheme 
with a monotonic flux limiter of the Osher's form. A Carte- 
sian coordinate system was used in the calculations. Their 
results showed that no hot spot was formed. Furthermore, 
disc formation was inhibited for higher 7 and no spiral 
sho ck was seen. 

iMakita. Mivawaki fc Matsudal il987f) carried out two- 
and three-dimensional numerical simulations of an accre- 
tion disc in a close binary system with a higher resolution 
using Simplified F lux vector Splitting finite volume method 
Jjvounouchi et alJll993t IShima fc Jvounouchilll99i) . Their 
computational region only covered the vicinity of the mass 
accreting star. The gas from the mass losing star was as- 
sumed to flow into the computational domain through a 
rectangular hole located at the LI Lagrangian point. They 
obtained quasi steady solutions, and confirmed the existence 
of both an accretion disc and spiral shocks for all cases with 
7 = 1.2, 1.1, 1.05 and 1.01. It was shown that a smaller 
7 value resulted in a more tightly wound spiral shock in 
two-dimensional calculations, but such tendency was not so 
obvious in three-dimensional cases. It was also shown that 
the spiral shock waves disappeared when tidal force due to 
the mass losing star was a rtificially cut off. 

iFuiiwara et alJ <l200lf) carried out three-dimensional 
simulations of a close binary system containing both the 
mass accreting star and the mass losing star, and investi- 
gated the interaction between the LI stream and the accre- 
tion disc. In the calculations, the ratio of constant specific 
heat was chosen as 7 = 1.2 and 1.01. They confirmed the 
existence of accretion disc as well as spiral shocks. No hot 
spot was found in the accretion disc. Instead, they found 
that a bow shock wave was formed due to the collision of 
the LI stream and the rotating disc. It was shown that this 
bow shock wave heated the outer part of the accretion disc, 
and also enhanced the density perturbation in the disc re- 
sulting in a more effective transfer of angular momentum by 
the tidal torque. 

From above numerical results, one can assume 
the existence of spiral shocks rather firmly even for 
three-dimensional accretion discs. This view is fur- 
ther supported by a r e cent observational evidence. 
ISteeehs. Harlaftis fc Hornel (^997) found the first convinc- 
ing evidence for spiral structure in the accretion disc of the 
eclipsing dwarf nova binary IP Pegasi using the technique 
known as Doppler tomography. 

Next, let us examine the question (ii). In the past nu- 
merical simulations of accretion discs, it has been customary 
to assume an adiabatic gas with a lower specific heat ratio 
than that of a mono-atomic gas (= 5/3) to account for non- 



adiabatic process of radiative cooling. So far, no simulation 
that accounts for radiative cooling effect with 7 = 5/3 has 
been made to show whether a disc structure is really formed 
and spiral shocks appear. It is also yet to know how good is 
the approximation to employ an adiabatic gas with a lower 7 
for simulating actual flowfield with radiative cooling effect. 

Recently. IStehlel l|l999f) carried out two-dimensional hy- 
drodynamic calculations of accretion discs including both 
the a-type viscosity and the effect of energy loss from the 
surface by radiation. They followed the ev olution of the lo- 
cal dis c thickness in the one-zone model of IStehle fc Spruitl 
(1999). The computational domain, however, only covered 
the vicinity of the mass accreting star. In their calculation, 
spiral shocks were clearly observed which dominated the disc 
evolution in a hydrodynamical time-scale. 

The purpose of this work is first to obtain the two- 
dimensional flowfield in a close binary system with radiative 
cooling effect, and then to explore whether a disc structure 
is really formed and spiral shocks appear in the accretion 
disc. Moreover, we would like to examine how good is the 
approximation of using a lower 7 value for simulating non- 
adiabatic processes of radiative cooling in the hydrodynamic 
calculation of accretion discs. We therefore solve the two- 
dimensional Euler equations with a constant specific heat 
ratio of mono-atomic gas (7 = 5/3) for a thin disc on the 
equatorial plane, and include radiative cooling effect by dis- 
charging energy in the vertical directions from the top and 
bottom surfaces of the disc. We employ the flux-limited dif- 
fusion (FLD) approximation to obtain radiative heat flux 
values. In the present study, we consider a realistic case of 
the observed close binary system. 

This paper is organized as follows. In ^5|we briefly de- 
scribe the numerical method, physical assumptions and nu- 
merical treatment of radiative transfer. In Sj^|we show the 
obtained results for accretion disc with radiative cooling. A 
detailed comparisons are made with the results of assum- 
ing an adiabatic gas with lower 7 values. In [jlj concluding 
remarks will be given. 



2 ASSUMPTIONS AND METHOD OF 
CALCULATIONS 

We consider a SU U Ma type of CV: OY Car as a realistic 
close binary system jRitter fc Kolblll995l) . The separation 
between the centers of these two stars is 4.22 x 10 10 [cm]. The 
orbital period is 1.55[hour]. The mass of these two stars are 
Mi = O.69M0 for the mass accreting star and M2 = O.O7M0 
for the mass losing star, respectively. In the present study, 
we assume the surface density of the mass losing star to be 
pa = 0.5 x 10 _9 [g cm" 3 ]. 



2.1 BASIC EQUATIONS AND PHYSICAL 
CONDITIONS 

The basic equations describing a two-dimensional flow of 
per fect gas in a rotating frame of refere nce can be described 
as iSawada. Matsuda fc Hachisul ll984') 



9Q + dF + dG +fl 
dt dx dy 



0, 



(1) 
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where p, u, v, e, p, f x , f y and F z are the density, the x and 
y components of the velocity, the total energy per unit vol- 
ume, the gas pressure, the x and y components of the net 
force that includes the gravity force, the Coriolis force and 
the centrifugal force, and the z component of the radiative 
energy flux, respectively. We ignore the radiation pressure 
because it is negligibly small if compared with the gas pres- 
sure. 

The equation of state is given by 



P : 



(7-l){e-f (u 2 +v 2 )} 



(3) 



which is that of an ideal gas characterized by a constant 
ratio of specific heats 7. 

A finite volume approach is used to discretize the gov- 
erning equations. The AUSM- DV scheme is used t o obtain 
the numerical convective flux l|Wada fc Lioulll994l) . In the 
time integration, a mat rix free LU-SGS method is employed 
jNakahashi et alJll999^ . Inner iterations are made to assure 
the s econd order temporal accu racy in the implicit integra- 
tion jYamamoto fc Kanol[l99rtl . A CFL number of 1000 is 
assumed in the calculations. 

The parameters that characterize the gas ejected from 
the mass losing component are the sound velocity, which is 
chosen to be Co = 73[km s -1 ], and the velocity of the gas 
inside of the star, which is uo = 0. Although we assume a 
stationary gas, the higher pressure inside of the mass losing 
star drives the gas flown into the computational domain. 
Note that if we employ AQ as a typical velocity scale as 
was in the previous works, where A is the separation of the 
two stars and Q is the angular velocity of the system, the 
dimensionless sound velocity of the ejected gas becomes 0.15. 
As a boundary condition at the surface of the mass accreting 
star as well as at the outer boundary of the computational 
domain, we assume a vacuum condition in which any gas 
reaching there is simply absorbed in. As an initial condition, 
we assume a gas with very low density, 0.5 x 10 _14 [g cm -3 ], 
and high temperature, 3.22x 10 6 /7[K], that fills up the entire 
computational domain. 

In order to avoid the numerical instability at the initial 
stage, we first calculate the flow of the gas with 7 = 5/3 
without radiative cooling up to 10 5 time steps using first- 
order-accurate scheme. After that we switch to second- 
order-accurate scheme. At the same time, the 7 value is 
changed for relevant cases. The radiative cooling effect is 
also included for non-adiabatic cases. 

Our computational domain includes both the mass los- 
ing star and the mass accreting star on the equatorial plane 
(Fig. . The computational domain that contains the mass 



losing star has 20 x 61 grid points, while that contains the 
mass accreting star has 44 x 61 grid points. These two do- 
mains are connected by a narrow strip that has 3 x 19 grid 
points. 



2.2 NUMERICAL TREATMENT OF 
RADIATIVE TRANSFER 

It is assumed that energy discharge due to radiative cooling 
occurs in the vertical directions from the top and bottom 
surfaces of the disc. Because the present flowfield is con- 
fined in a two-dimensional equatorial plan, it is necessary to 
assume the distribution of physical quantities in the vertical 
direction in the disc for calculation of radiation. We assume 
a hydrostatic balance in the vertical direction to determine 
the thickness of the disc, and that the physical quantities in 
the vertical direction are constant with the local values in 
the equatorial plan. The thickness of the disc H is given by 



H 



K 



GM\ 



(4) 



where c s , f2j^, n and G are the sound speed of the gas on 
the disc, the Keplerian frequency, the radial distance from 
the mass accreting star and the gravitational constant, re- 
spectively. 

The coupling of the flowfield and the radiation field is 
occurred through the term of dF z /dz that appears in the 
energy equation. dF z /dz is obtained by solving the following 
equation for radiation energy 

= 4ir X B - c X E, (5) 

where B, E, c and x are the Planck function, the radiation 
energy density per unit volume, the speed of light and opac- 
ity, respectively. Because the speed of light is far larger than 
the sound velocity, and we only seek for a steady state so- 
lution, we ignore the time differential term that appears in 
the radiative transfer equation. 

In the present work we employ the Kramer's opac- 
ity which has no dependence on frequency. This Kramer's 
opacity is suitable for the calculation of pure hydrogen 
;as above 10 4 [K1, and is written in the follo wing form 
Cannizzo. Shaffer fc Wheeledll98gUStehldll999l) 

X = 2.8x 10 2 V 2 T -35 [cm -1 ], (6) 

where T is the temperature of the gas. 

As mentioned earlier, we consider the radiative transfer 
only in the vertical direction of the disc, i.e. the z direction 
perpendicular to the x-y plane. The radiative transfer within 
the disc is ignored. Moreover, the radiative transfer is taken 
into account in the calculations only in the computational 
domain that contains the mass accreting star and hence the 
accretion disc. 

The radiative heat flux in the accretion disc is cal- 
culated by the FLD approxim ation iAlme fc Wilson! Il974l : 
iLevermore fc Pomraning|l98ll) . We assume the gas is in lo- 
cal thermodynamic equilibrium at a temperature that needs 
not correspond to that of the radiation field. This is because 
we solve the radiation energy equation to determine the 
radiation energy density E. With the FLD approximation, 
the radiativ e flux can be written in the fo rm of Fick's law 
of diffusion jLevermore fc Pomranin5 ll9811 as 
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Figure 1. The generalized curvilinear coordinate used here. The mesh size is 20 X 61 in the mass losing star region, 3 X 19 in the 
connecting region and 44 X 61 in the mass accreting star region. 



-D 



dE_ 

dz ' 



with a diffusion coefficient, D, given by 

X 



(7) 



(8) 



The dimensionless function A = X(E) is called as th e flux 
limiter. We employ Minerbo's model (Mincrbo 1978) as a 
choice of A, i.e. a constraint on the anisotropy of the radi- 
ation field. Minerbo assumed a piecewise linear variation of 
the specific intensity with angle, and found the functional 
form of the flux limiter as 



A(i?) 



2/(3 + V9 + 12R 2 ) 
(l + R+Vl + ZRy 1 



< R «S 3/2, 
3/2 < R < oo, 



(9) 



where R is a dimensionless quantity R = \dE/dz\ /{%E). In 
the optically thin limit R — > oo, the flux limiter gives 



lim \(R) = -, 

-R— ►oo i\, 



(10) 



to the first order in R~ . The magnitude of the flux therefore 
approaches \F Z \ = c\dE/dz\ /(xR) — C E, which obeys the 
causality constraint. In the optically thick or diffusion limit 
R — * 0, the flux limiter gives 

1 

3' 



lim \(R) 

R->0 



(11) 



to the first order in R, so that the flux takes the value given 
by equation J7J . 

The method of the solving the radiation energy equation 
JSJ is described in Appendix 1X1 



3 RESULTS 

Figs. |U (a), (b) and (c) show the density distribution in a 
(quasi) steady state (at 10 6 step) for the cases of adiabatic 
gas with 7 = 5/3, 1.2 and 1.01, respectively. In these calcu- 
lations radiative cooling effect is not included. It is shown 
that the disc structure is formed and the spiral shocks are 
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Figure 2. Density contours for each case, (a), (b) and (c) show the density contours for the cases of 7 = 5/3, 1.2 and 1.01 not accounting 
for the radiative cooling, respectively, (d) shows the density contours for the case of 7 = 5/3 accounting for the radiative cooling. The 
density range in (a), (b), (c) and (d) is from 1.0 X 10~ 13 to 1.0 X 10~ 12 , from 1.7 X 10~ 13 to 2.6 X 10~ 12 , from 5.4 X 10~ 13 to 1.1 X 10 -11 
and from 2.5 X 10~ 13 to 4.9 X 10 — 12 [g cm -3 ], respectively. Each density range is divided by 20 lines with an equal increment. The system 
rotates counterclockwise. 



appeared in each case. For the case of 7 = 5/3, the spiral 
arms are less clear as shown in Fig.|2|(a). This is because the 
contour lines are plotted at equal interval while the density 
in the disc for this case is relatively lower than that in the 
LI stream. 

It should be noted that the steady state is not steady in 
a strict sense, but slightly oscillates. Especially, for the case 
of 7 = 5/3, the density pattern and the shock position peri- 
odically change with time, although the general pattern in 
the inner and outer region remains unchanged. The temporal 
fluctuation of the disc structure, however, is not significant 
and the density distribution shown in Fig. [5] (a) gives a typ- 
ical one for the case of 7 = 5/3. This is suggested by the 
fact that the present density distribution virtually coincides 
with that of the ensemble average taken from 2.0 X 10 5 to 
1.0 x 10 6 time steps. 



The pitch angle of spiral arms has a clear cor- 
relation with 7. The lower 7 thus leads to a tightly 
wound spiral shocks. These results are consistent with 
the previous works JSaw ada. M atsuda fc Hachisul Il986al : 
iMakita. Mivawaki fc Matsudalll987l) . 

The maximum Mach number in the disc is about 4, 5 
and 23 for the cases of 7 = 5/3, 1.2 and 1.01, respectively. 

Fig. [5] (d) shows the density distribution for the case in 
which the specific heat ratio is fixed as 5/3 and the effect of 
radiative cooing is accounted for. Obviously, one can see that 
the disc structure is formed and spiral shocks are appeared. 

Spiral shocks disappear near the mass accreting star 
and nearly axisymmetric disc structure is formed. Those flow 
patterns are quite similar to that obtained in the calculation 
of adiabatic case with 7 = 1.01, although the absolute value 
of density in the disc becomes smaller. 
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Figure 3. Temperature contours with logarithmic scale for each case, (a), (b) and (c) show the temperature contours for the cases of 
7 = 5/3, 1.2 and 1.01 not accounting for the radiative cooling, respectively, (d) shows the temperature contours for the case of 7 = 5/3 
accounting for the radiative cooling. The temperature range in (a), (b), (c) and (d) is from log 10 T[K] = 5.9 to 8.1, from 5.9 to 7.5, from 
4.6 to 7.2 and from 4.6 to 7.2, respectively. The temperature range is divided by 20 lines with an equal increment. 



The maximum Mach number in the disc is found to be 
23 for this case which is the same as that for the case of 
7 = 1.01. 

Figs. |21 (a), (b) and (c) show the temperature distribu- 
tion in a (quasi) steady state with logarithmic scale for the 
cases of adiabatic gas with 7 = 5/3, 1.2 and 1.01, respec- 
tively. As can be seen in Fig. 0(a) and (b), the temperature 
of the gas in the disc is very high and becomes above 10 7 [K] 
near the mass accreting star. The smaller 7 value for this 
case results in nearly isothermal distribution in the disc that 
gives substantially lower temperature near the mass accret- 
ing star. As shown in Fig. (c) , compared with the other 
adiabatic cases, the temperature distribution for the case of 
7 = 1.01 is found to be close to flat. 

Fig. 0(d) shows the temperature distribution with loga- 
rithmic scale in a (quasi) steady state for the case of 7 = 5/3 
with radiative cooling. Compared with the other adiabatic 



cases without radiative cooling, the temperature in the ac- 
cretion disc is significantly lowered to less than 10 5 [K], ex- 
cept for the small disc region near the mass accreting star 
where the temperature exceeds 10 [K]. In other words, the 
accretion disc exhibits a dual structure, i.e. the significantly 
cooled outer region and the inner core region where temper- 
ature is very high. The precise structure in the inner core 
region, however, is difficult to resolve with the present mesh 
system, because the pitch angle of the spiral shocks becomes 
smaller as approaches to the inner region. 

The evolution of the total mass in the mass accreting 
region is shown for each case in Fig. The total mass is 
defined as a sum of all the masses in the mass accreting 
region. The horizontal axis denotes the time step. 1.0 x 10 6 
time steps correspond to more than 700 orbital periods. In 
all the cases, the total mass in the mass accreting region 
becomes almost constant after 3.0 x 10 time steps but begin 
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Figure 4. The evolution of the total mass in the mass accreting region denoted by solar mass unit. The total mass is defined as a sum 
of all the masses in the mass accreting region. 



to oscillate slightly around each mean value. Therefore, we 
can say that each flowfield in the accretion disc reaches a 
(quasi) steady state. The total mass for the case accounting 
for radiative cooling is found to be larger than that for the 
case of 7 = 1.2, and smaller than that for the case of 7 = 
1.01. From the results, we can say that the overall behavior 
of the flowfield as well as the evolution of the total mass of 
the adiabatic case with a smaller 7 value is quite similar to 
the one with radiative cooling shown in the present study. 
Although the precise 7 value can change for case-to-case, one 
can conclude that the use of a lower 7 value for simulating 
radiative cooling effect is almost a valid assumption. 

Finally, let us examine the influence of gas properties 
from the mass losing star. Because the process of discharge 
and absorption of radiation depends on density and tem- 
perature of the gas in the accretion disc that is influenced 
primarily by the gas from the mass losing star, we need 
to know what are the consequences of changing proper- 
ties of the gas. Fig. [3 (a) shows the density contours in 



the accretion disc for the lower temperature case in which 
po = 0.5 x 10~ 9 [g cm" 3 ] and Co = 48[km s" 1 ]. In the calcu- 
lation, the radiative cooling effect is considered. In this case, 
the accretion disc slightly shrinks, but the overall feature is 
unchanged if compared with the standard case shown in Fig. 
IEl(d). On the other hand, if we increase the surface density 
of the mass losing star, some changes are seen in the re- 
sults. Fig. 13(b) shows the density contours for the case with 
po = 0.5 x 10" 8 [g cm" 3 ] and c = 73[km s" 1 ]. Note that the 
density scale is changed to show the overall distribution. In 
this case, the inner disc region again appears where density 
becomes very low, but its radius becomes smaller. Surround- 
ing this inner region, there appear several spiral arms in the 
intermediate region. This intermediate region has a ring like 
structure, which is connected to an outer spiral arm and also 
to the incoming flow from the Llpoint. In Fig. [3(c), the den- 
sity contours for the limiting case of po = 1.0 x 10~ 7 [g cm" 3 ] 
and Co = 24[km s" 1 ] are shown. One can see that the overall 
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(a) : (density), p =0.5x1 0" a [g cm" 3 ], c =4S[km a" 1 ] (b) : (density), p =0.5x1 0^[g cm" 3 ], c fl =73 [km s 1 ] 
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Figure 5. Density contours for the case of 7 = 5/3 accounting for the radiative cooling for the different parameters of po and cq. 
(a), (b), (c) and (d) show the density contour for the cases of (po,c ) = (0.5 X 10 -9 , 48), (0.5 X 10 -8 , 73), (1.0 X 10 _7 ,24) and (0.5 X 
10 — 9 [g cm - 3 ],73[km s^ 1 ]), respectively. The figure of (d) is same as that of Fig.[2](d), which is the standard case. The density range 
in (a), (b) and (c) is from 2.5 X 10~ 13 to 4.9 X 10~ 12 , from 2.4 X 10~ 12 to 4.4 X 10~ n and from 4.7 X 10~ 12 to 8.4 X 10 _11 [g cm" 3 ], 
respectively. The density range in (d) is same as it of (b). Each density range is divided by 20 lines with an equal increment. 



feature is quite similar to that shown in Fig. [3(b), though 
the inner disc region further shrinks. 



4 CONCLUDING REMARKS 

Two-dimensional hydrodynamic calculations of an invisicid 
flowfield in a close binary system are carried out by solving 
the Euler equations with radiative transfer. In the present 
study, the specific heat ratio is assumed to be constant while 
radiative cooling effect is included as a non-adiabatic pro- 
cess. The cooling effect of the disc is considered by discharg- 
ing energy in the vertical directions from the top and bot- 
tom surfaces of the disc. We use the flux-limited diffusion 
approximation to calculate the radiative heat flux values. 
The obtained results are summarized as follows: 



(i) Around the mass accreting star, there appears an ac- 
cretion disc and spiral shocks are formed in the disc, even 
for the case with radiative cooling. 

(ii) The total mass of the disc reaches a (quasi) steady 
value. It becomes larger when radiative cooling effect is ac- 
counted for than that assuming 7 = 5/3, and is close to the 
value for the case of lower 7. 

(iii) The lower 7 value for an adiabatic gas which has been 
employed in the simulation of accretion discs in order to 
account for radiative cooling effect is found to be almost a 
valid assumption in terms of the overall flow features and 
also the amount of the total mass. 

What are shown in this work are the results of two- 
dimensional hydrodynamic calculation of inviscid flowfield 
accounting for radiative cooling effect. In order to obtain 
more rigorous results for the accreting flowfield, we need 
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to explore three-dimensional flowfield that is coupled with 
radiative transfer. In such case, one needs to account for 
radiative transfer that is parallel to equatorial plane. This 
is because the accretion disc is not necessarily a flat thin 
disc. It is also needed to employ a more detailed radiation 
model that can consider spectral dependence of opacity. For 
such calculation, a flux-limited diffusion model is no more 
applicable and we need further to employ a detailed simu- 
lation method such as to use ray-tracing technique to solve 
radiative transfer equations. Our future studies will focus on 
these aspects. 
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APPENDIX A: DISCRETIZING OF 
RADIATION ENERGY EQUATION 

Let us consider the method of solving the radiation energy 
equation JSJ . Using the flux- limited diffusion approximation 
10, equation © can be rewritten in the following form for 
radiation energy density E as 

Because the distribution of physical quantities in the vertical 
direction is assumed to be constant, the diffusion coefficient 
D is independent of the height z. Therefore, we replace the 
partial differential in z by multiplying 1/H. As a result, 
equation HA 111 is discretized as 

2 ° n jWyI = ^* Rn ~ cx n E n+1 , (A2) 

where the superscript n denotes the time step. The factor 2 
in the left hand side of the equation (IA2II appears because 
discharging energy by radiation occurs from the top and bot- 
tom surfaces of the disc. Except for the radiation energy den- 
sity E, the dependent variables are evaluated at time step 
n. Equation JX2J is solved for E n+1 . This E n+1 is substi- 
tuted into the discrete form of dF z /dz, i.e. 2D n E n+1 /(H n f 
which is further substituted into the energy equation of the 
Euler equations. This completes the coupling procedure of 
the flowfield with radiation. 

This paper has been typeset from a T^rX/ DT^rX file prepared 
by the author. 



